In [ ]:
import os
import sys
# Third-party
from astropy.io import ascii
import astropy.coordinates as coord
import astropy.table as table
import astropy.units as u
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as pl
pl.style.use('apw-notebook')
%matplotlib inline
from astropy.io.fits import getdata
import gary.coordinates as gc
import gary.dynamics as gd
import gary.potential as gp
from gary.dynamics.orbit import combine as combine_orbit
from gary.dynamics.core import combine as combine_ps
from gary.observation import distance_modulus
from ophiuchus import galactocentric_frame, vcirc, vlsr
import ophiuchus.potential as op
In [ ]:
catalog_data_path = "/Users/adrian/projects/globber/data/gc_catalogs/"
cluster_name = 'NGC 5897'
nsamples = 128
potential = op.load_potential('static_mw')
In [ ]:
pm_gc_main = np.genfromtxt(os.path.join(catalog_data_path,"gl_2012_J2000.cat1.txt"), dtype=None,
skip_header=2,
usecols=[0,2,3,6,7,8,9,10,11,12,13],
names=['ngc_num','ra','dec','dist','dist_err','mu_ra','mu_ra_err',
'mu_dec','mu_dec_err', 'vr', 'vr_err'])
pm_gc_main = table.Table(pm_gc_main)
go = ascii.read(os.path.join(catalog_data_path,"go97_table1.txt"))
all_gc = pm_gc_main
all_gc['name'] = np.array(["NGC {}".format(x) for x in all_gc['ngc_num']])
all_gc = table.join(all_gc, go, keys='name')
# HACK: use distance modulus from Kock et al. and make uncertainty better
from gary.observation import distance as distance_from_dm
cluster = all_gc[all_gc['name'] == cluster_name]
cluster['dist'] = distance_from_dm(15.55).to(u.kpc).value
cluster['dist_err'] = 0.05*cluster['dist']
cluster
In [ ]:
cluster_c = coord.ICRS(ra=cluster['ra']*u.degree,
dec=cluster['dec']*u.degree,
distance=cluster['dist']*u.kpc)
In [ ]:
xyz = cluster_c.transform_to(galactocentric_frame).cartesian.xyz
vxyz = gc.vhel_to_gal(cluster_c, pm=(cluster['mu_ra']*u.mas/u.yr,
cluster['mu_dec']*u.mas/u.yr),
rv=cluster['vr']*u.km/u.s,
galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
In [ ]:
np.random.seed(42)
_distances = np.random.normal(cluster['dist'], cluster['dist_err'], size=nsamples)
cluster_samples_c = coord.ICRS(ra=(np.zeros(nsamples) + cluster['ra'])*u.degree,
dec=(np.zeros(nsamples) + cluster['dec'])*u.degree,
distance=_distances*u.kpc)
_mu_ras = np.random.normal(cluster['mu_ra'], cluster['mu_ra_err'], size=nsamples)
_mu_decs = np.random.normal(cluster['mu_dec'], cluster['mu_dec_err'], size=nsamples)
_vrs = np.random.normal(cluster['vr'], cluster['vr_err'], size=nsamples)
# ---
samples_xyz = cluster_samples_c.transform_to(galactocentric_frame).cartesian.xyz
samples_vxyz = gc.vhel_to_gal(cluster_samples_c,
pm=(_mu_ras*u.mas/u.yr, _mu_decs*u.mas/u.yr),
rv=_vrs*u.km/u.s,
galactocentric_frame=galactocentric_frame,
vcirc=vcirc, vlsr=vlsr)
In [ ]:
w0 = gd.CartesianPhaseSpacePosition(pos=xyz, vel=vxyz)
mean_orbit = potential.integrate_orbit(w0, dt=-0.5, nsteps=12000)
_w0 = gd.CartesianPhaseSpacePosition(pos=samples_xyz, vel=samples_vxyz)
all_w0 = combine_ps((w0,_w0))
orbit = potential.integrate_orbit(all_w0, dt=-0.5, nsteps=12000)
In [ ]:
w0.represent_as(coord.CylindricalRepresentation)
In [ ]:
c_back,v_back = orbit[:100].to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit -50 Myr
galactocentric_frame=galactocentric_frame)
orbit_forw = potential.integrate_orbit(all_w0, dt=0.5, nsteps=100)
c_forw,v_forw = orbit_forw.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr, # orbit +50 Myr
galactocentric_frame=galactocentric_frame)
In [ ]:
fig,axes = pl.subplots(8,8,figsize=(12,12),sharex=True,sharey=True)
for i in range(len(axes.flat)):
ax = axes.flat[i]
# put markers at tidal radius along radial vector from center of MW
this_xyz = samples_xyz[:,i]
r = np.linalg.norm(this_xyz,axis=0)
r_hat = (this_xyz / r)
r_tide = (cluster['M'][0] / potential.mass_enclosed(this_xyz))**(1/3.) * r
pos,neg = this_xyz + r_hat*r_tide, this_xyz - r_hat*r_tide
Lpts = coord.Galactocentric(coord.CartesianRepresentation(np.vstack((pos.value,neg.value)).T*u.kpc))
Lpts_icrs = Lpts.transform_to(coord.ICRS)
ax.set_xticklabels([])
ax.set_yticklabels([])
_style = dict(ls='-', alpha=0.5, lw=2., marker=None)
ax.plot(c_back.ra.degree[:,i], c_back.dec.degree[:,i], color='k', **_style)
ax.plot(c_forw.ra.degree[:,i], c_forw.dec.degree[:,i], color='#2166AC', **_style)
ax.plot(cluster_c.ra.degree, cluster_c.dec.degree, marker='o', color='k')
ax.plot(Lpts_icrs.ra.degree, Lpts_icrs.dec.degree, marker='o', color='r', ls='none')
ax.set_xlim(235,225)
ax.set_ylim(-25,-15)
fig.tight_layout()
In [ ]:
cluster
In [ ]:
from gary.dynamics.mockstream import fardal_stream, streakline_stream
import gary.integrate as gi
In [ ]:
stream_models = dict()
In [ ]:
fig,axes = pl.subplots(4,4,figsize=(12,12),sharex=True,sharey=True)
for i in range(len(axes.flat)):
print(i)
ax = axes.flat[i]
if i in stream_models:
stream = stream_models[i]
else:
prog_orbit = potential.integrate_orbit(all_w0[i], dt=-0.25, nsteps=2000)
stream = fardal_stream(potential, prog_orbit[::-1], prog_mass=cluster['M'][0],
release_every=1, Integrator=gi.DOPRI853Integrator)
# stream = streakline_stream(potential, prog_orbit[::-1], prog_mass=cluster['M'][0],
# release_every=2., Integrator=gi.DOPRI853Integrator)
stream_models[i] = stream
stream_c,stream_v = stream.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr,
galactocentric_frame=galactocentric_frame)
ax.plot(stream_c.ra.degree, stream_c.dec.degree, marker='.', alpha=0.1, ls='none')
ax.set_xlim(235,225)
ax.set_ylim(-25,-15)
fig.tight_layout()
In [ ]:
fig,axes = pl.subplots(4,4,figsize=(12,12),sharex=True,sharey=True)
for i in range(len(axes.flat)):
ax = axes.flat[i]
stream = stream_models[i]
stream_c,stream_v = stream.to_frame(coord.ICRS, vcirc=vcirc, vlsr=vlsr,
galactocentric_frame=galactocentric_frame)
ax.plot(stream_c.ra.degree, distance_modulus(stream_c.distance), marker='.', alpha=0.1, ls='none')
ax.set_xlim(235,225)
ax.set_ylim(15,16)
fig.tight_layout()
In [ ]: